      subroutine jet(theta, rminor) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      real(double) , intent(in) :: theta 
      real(double) , intent(inout) :: rminor 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      real(double) :: rmaj, s, a, b, c, sint, cost, alpha, beta, gamma, delta1&
         , delta2, delta3, gl, dr, counter, glp 
!-----------------------------------------------
!
      data rmaj/ 2.96/  
      data s/ 1./  
      data a/ 0.119/  
      data b/ 0.0579/  
      data c/  - 0.3064/  
!
!     given theta, calculate rminor by bisection
!
      sint = sin(theta) 
      cost = cos(theta) 
      alpha = a*sint**2 + b*cost**2 
      beta = 2.*rmaj*a*sint 
      gamma = a*rmaj**2 
      delta1 = c*sint**2 
      delta2 = 2.*c*rmaj*sint 
      delta3 = c*rmaj**2 
!
      gl = (alpha*rminor**2 + beta*rminor + gamma)**2 + delta1*rminor**2 + &
         delta2*rminor + delta3 + s 
!
      dr = 0.1*rmaj 
      counter = 0.0 
      rminor = 0. 
!
      rminor = rminor + dr 
      if (rminor < 0.) then 
         rminor = 0. 
         dr = -dr 
      endif 
      counter = counter + 1. 
!
      glp = (alpha*rminor**2 + beta*rminor + gamma)**2 + delta1*rminor**2 + &
         delta2*rminor + delta3 + s 
!
      if (glp*gl < 0.0) dr = -0.5*dr 
      gl = glp 
      do while(abs(gl)>1.E-5*rmaj .and. counter<100.) 
         rminor = rminor + dr 
         if (rminor < 0.) then 
            rminor = 0. 
            dr = -dr 
         endif 
         counter = counter + 1. 
!
         glp = (alpha*rminor**2 + beta*rminor + gamma)**2 + delta1*rminor**2 + &
            delta2*rminor + delta3 + s 
!
         if (glp*gl < 0.0) dr = -0.5*dr 
         gl = glp 
      end do 
!
      return  
      end subroutine jet 
